function [output,T1,T2] = fn_tax(delta,R1,D0,chi,s,phi,dist_d,d_min,d_max)

rho1_s = fn_rho1(s,phi);
chi_s  = fn_chi(s,chi);

f = @(x) min(x*R1,delta).*pdf(dist_d,x);

T1 = integral(f,d_min,d_max);
T2 = chi_s.*rho1_s*D0;
T  = T1-T2;

output = max(T,0);

end